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Abstract 

Populations of uncoupled limit-cycle oscillators receiving common random impulses show various 
types of phase-coherent states, which are characterized by the distribution of phase differences 
between pairs of oscillators. We develop a theory to predict the stationary distribution of pairwise 
phase difference from the phase response curve, which quantitatively encapsulates the oscillator 
dynamics, via averaging of the Frobenius- Perron equation describing the impulse-driven oscillators. 
The validity of our theory is confirmed by direct numerical simulations using the FitzHugh-Nagumo 
neural oscillator receiving common Poisson impulses as an example. 
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I. INTRODUCTION 
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Coherence phenomena exhibited by dynamica 
been the focus of much recent research l|, [3, ^ 
Experimentally, synchronization among dynamical units receiving common fluctuating drive, 
or response reproducibility of a single unit receiving identical fluctuating drive, has been 
shown in neurons 1, 2], 3|, chaotic lasers 4], and electrical oscillators 0, 1?]. The slightly 
counterintuitive phenomenon of desynchronization or anti-reliability via a common input 
has been seen in electrical oscillators Q], electrochemical oscillators [8], and light-sensitive 
circadian cells j^. Further, coexistence of multiple synchronized groups of dynamical units 
have been observed in chaotic electrical circuits, and are known as multiple basins of consis- 
tency [6| . For limit-cycle oscillators, theoretical analysis has yielded quite a few quantitative 
results explaining synchronization, desynchronization, and multiple syrichronize d g roups or 
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clusters exhibited in an ensemble of limit-cycle oscillators 

m 

Our previous work [3, [isl analyzed the linear stability of synchronized or clustered states 
of uncoupled limit-cycle oscillators subject to random common external impulses by calcu- 
lating the Lyapunov exponent, which quantifies the average rate of growth of an infinitesimal 
phase separation between a pair of oscillators. The only dynamical information we require 
about the oscillator is contained in a simple function called the phase response curve (PRC) 
describing the magnitude of phase advance or retardation due to a perturbation at a given 
phase jiol . [2^. The PRC has been measured in many oscillator- like systerns,_ including 



neurons, circadian oscillators, cardiac cells, and electrical circuits 
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non-frequent impulses, the Lyapunov exponent A is given by 



A = A 





b J dc\n 


/o 





l + -^G{<f),c) 



p{c) 



2^. For 
(1) 



where A is the mean number of impulses in a unit time (or rate), Gi(j),c) is the PRC for 
an impulsive perturbation whose intensity and direction (or mark [25]) is c, p{c) is the 
probability density of the mark, and the integral is over the oscillator phase and the mark 
c. A negative (positive) A means that an infinitesimal phase difference shrinks (grows) on 
the average, resulting in synchronization (desynchronization) of the oscillators. 

However, the Lyapunov exponent alone is not sufficient in characterizing the whole coher- 
ence phenomena induced by the common impulses, because it is an average quantity over the 
entire limit cycle that characterizes only the local linear stability of the synchronized state. 



The phase difference generally does not monotonically decrease or increase over successive 
common impulses due to ffuctuations in the expansion rates of the phase difference, which is 
determined by the precise form of the PRCs. When small external noises or inhomogeneities 
exist, such fluctuations may induce large excursions from the synchronized state even if the 
Lyapunov exponent is negative on average. Oscillator pairs may find themselves with large 
phase difference, but the global distribution of the phase difference cannot be explained by 
a linear stability analysis. 

In this paper, we further the theoretical analysis for an ensemble of generic uncoupled 
limit-cycle oscillators to obtain the stationary distribution of pair- wise phase difference 36 1. 
Starting from general dynamical equations for a pair of limit-cycle oscillators driven by com- 
mon impulses, we derive a pair of random maps and the corresponding two-body Frobenius- 



Perron equation 
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27l | using the phase reduction method [7|, ll9|, l20|. We then derive an 



approximate one-body Frobenius-Perron equation for the phase difference by averaging out 
the fast phase dynamics, which yields the stationary distribution of the phase difference. The 
theoretical result is compared with direct numerical simulations using FitzHugh-Nagumo os- 
cillators receiving common Poisson impulses. 



II. THEORY 



A. Phase reduction of the dynamical equation 

We investigate a pair of uncoupled oscillators receiving common random impulses and also 
subject to a weak additive Gaussian white noise independently. The stochastic dynamical 
equation for the i-th oscillator in this pair is 

N{t) 

X,{t) = F{X,) + c^^^)h{t - + y/DH{Xi)rj,, (2) 

ra=l 

where i = 1,2, Xi{t) G R'^^ is the oscillator state at time t, F{Xi) : R'^' R^' the 
dynamics of a single oscillator, N(t) the number of received impulses up to time t, t*^"-' the 
arrival time of the n-th impulse, c*-"-* G R^ the intensity and direction, or mark [25|, of the 
n-th impulse, (T{Xi, c) : R^ x R^ —>■ R^'^ is the coupling function describing the effect of 
an impulse c to Xi, h(t — t^"'^) is the infinitesimally narrow unit impulse whose waveform 
is localized at the time t^") of the impulse (/^ - t^''^)dt = 1), H{X,) G i?^^x^^ the 
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coupling matrix of the independent noise to the oscillator, rji G -R^^ a Gaussian white 
noise of unit intensity with correlation {f]fit)Vj{s)) = ^(^ ~ s)6ai3Sij added independently to 
each oscillator, and D the intensity of the independent noise. We interpret Eq. ([2]) in the 
Stratonovich sense. If the impulses and the independent noises are absent {H = 0, cr = 0), 
the system is assumed to have a single stable limit-cycle solution, Xo(t). 



As in our previous papers 0, ll6|, we use the phase reduction method to 



analyze the 



dynamics of impulse-driven oscillators. We define an asymptotic phase [19|, |20| </> along 
the limit cycle Xo{t) that constantly increases with a natural frequency u, and extend the 
definition of phase to the whole state space of the oscillator (except phase singular sets) 
by identifying the orbits that asymptotically converge to the same point on the limit cycle. 
This defines a mapping from the oscillator state X G -R^^ to the phase (p ^ [0? !]• 

We assume that the interval between impulses is long compared to the relaxation time 
back to the limit-cycle, so the oscillator is almost always on the limit-cycle when an impulse 
is received. We can then reduce Eq. (|2]) to the dynamics of a single asymptotic phase (pi. The 

(n) 

dynamics of the phase (f)\ right before the n-th impulse is received can be approximately 
described by a random map 

where G'(0, c) is the PRC, tur*^"^ is the increase in phase during the interval between the 
n-th and (n + l)-th impulses r'-"^ = — t^"-), and 7^"^ is the displacement caused by the 

additive independent Gaussian noise rji in the interval t^"-\ From now on, we assume that 
the range of (p to be the real numbers R by taking into account the number of windings 
around the limit cycle, which makes the treatment of periodic boundary conditions easier 



in the following derivation 28 1 



The PRC G{(p, c) describes the change in phase of the oscillator when an impulse of mark 
c is received at phase (p on the limit cycle, which is periodic in (p, i.e. G(6+ 1, c) = G{(p, c). 



It can be obtained 
term in Eq. ([2]) as 



jy applying the approximation theorem by Marcus 2l(| to the impulsive 
7| 



c) = (Xo(0) + g{Xo{<P), c)) - 0, (4) 

I — I 

where g{X,c) = jexp (y^j{X,c){d/dXj)^ - ijx [37]. The PRC is related to the 



phase sensitivity function [20] Zi{(p) = c^0/f?-^ilx=Xo(</>) ^(0) — -^(0) ' '^(^o(0))C) 
when the effect of the impulse (t(Xo(0), c) is small. 



Generally speaking, the displacement 7^"^ depends on the oscillator phase ^j-"'^ the im- 
pulse mark c^"\ and the relaxation path to the limit cycle after each impulse. We approx- 
imate the actual distribution function of 7^^"'' by a zero-mean Gaussian normal distribution 
with variance e'^r^^^ [38i]. The approximate diffusion constant e can be obtained by ignoring 
the fast relaxation dynamics to the limit cycle after the impulse and by averaging the phase 



dependence over the limit cycle as 



= ['Y1 U<P)Zj{<P)H{Xo{<P))ikH{Xo{<l)))jkd<l), (5) 

where we utilize the fact that the stationary phase distribution of a single oscillator receiving 
infrequent impulsive forcing is nearly uniform . As we demonstrate later, this is a good 

approximation for oscillators whose relaxation to the limit cycle is sufficiently fast. 



B. Frobenius-Perron equation for the phase difference 

Let us consider the dynamics of the joint probability distribution p(0i, 02, n) of the phases 
(015 02) right before the n-th impulse, determined by the random map Eq. ([3]). We assume 
the range of phase variables to be 0i,2 ^ R- The Frobenius-Perron equation for the evolution 
of the joint distribution is 

p(0i,02,n + 1) 

/OO /"OO POO p noo POO 

d(P[ d<P'^ dr dc d-ii rf72Vr(r)p(c)i?(7i,r)i?(72,r) X 
-00 J —00 Jo Jc J ~oo J ~oo 

5 (01 - 01 - ^(01, C)-UJT - 7l)5(02 - 02 - ^(02, c) - UJT - I2) p{,<P'l, 02, n) 

/OO POO POO P 

d(i)[ / rf02 dr dcW{T)p{c)R{(j)i - - G{(f)[, c) - lut, r) x 
-00 J —00 Jo J c 

^(02 - 02 - ^(02, C) - UJT, T)p{(f)[, 02, n), (6) 

where W{t) is the inter-impulse distribution, G(0, c) is the PRC, and -R(7i, r) is the probabil- 
ity that an oscillator i has diffused an amount 7^ in a time interval r, which we approximated 
as a normal distribution with variance e^r. 
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Going to the center-of-mass coordinates, we change variables to ip = (0i + 02)/2 and 
= (pi —02, where ip is the mean phase and ^ is the phase difference. The Frobenius-Perron 
equation is transformed as 

/oo poo roo p 

di)' / di' dr dcp{c)W{r)x 
-oo J — oo ^0 J c 

Riip + l-ip' -^-Ciip' + ^,c] -ujt,t] X 



2 2 y 2 ^ 

- f - ^' + I - G (^^' - |,c^ - LUT,r^ p{^',e,n). 

We now restrict the mean phase to G [0, 1) and the phase difference to ^ G (—1, 1) 
similarly to Ermentrout and Saunders |28|] by introducing a new distribution function 



oo oo 



P{^P,^,n)=Y^ 5^ p(^+P,e + 2g,n), (7) 

p=— oo (?=— oo 

which sums up contributions from pairs of phase values with different winding numbers 
but represent physically equivalent situations on the limit cycle. This "wrapped" P{ip,^,n) 
corresponds to the actual distribution of the mean phase and the phase difference measured 
in simulations or experiments. Using the periodicity of the PRC, we obtain 

/•l /*1 pco p 

P(^,e,n+1)= d^p' di' i dr dcp{c)W{T)x 

R(^^-i-^> + ^ + q_G(^^P'-^,c^ - cor,r^ P{^P', e', n) , 

where the summation involves all pairs of p and q of equal parity (vr(-) denotes the parity 
of an integer). 

To obtain a closed equation for the phase difference ^, we now average out the fast 
dynamics of the mean phase, ip. If the impulses are not so frequent and the magnitude of 
the independent noise is small, the mean phase ^ is a rapidly changing variable compared 
to the phase difference ^. Then ip and ^ can be taken to be nearly independent, and the 
joint probability density can be separated as P{ip, ^, n) ~ S{ip, n)U{^, n), where S{ip, n) and 
U (^, n) are the probability density functions of ip and ^, respectively. Note that U (^, n) 
is periodic in ^, U{S, ± l,n) = U{^,n), because ^ and ^ ± 1 represent the same phase 



difference. For non-frequent impulses, ip is almost uniformly distributed on the limit cycle, 
S{t/j, ^ 1 0, [isl. We then average over the tjj on both sides to obtain 



U{^,n + 1) = I d^' I dij' I 
J-i Jo Jo 



d^T{^,i,^lj\OU{i\n), (8) 



where 



dr j dcp{c)W{r)yi 



^ ( ^ - I - V^' + I + g - G f ^' - |, c ) - a;r, r ) . (9) 



We now derive an approximate one-body Frobenius-Perron equation for the distribution of 
the phase difference 

U{i, n+l) = j^^ OU{i\ n)d^', (10) 
where the transition probability is given by 

^(^,0 = r#' r#T(^,e,^',o- (11) 

Jo Jo 

Namely, we have reduced the problem to finding the stationary distribution of a Markov 
process for the random variable ^ with transition probability X{C,,C,')- By numerically es- 
timating the transition probability X(^,^') from the PRC, Eq. (fTOj) can be iterated until a 
stationary state is reached. X{^,^') is periodic in ^ and ± 1,^' ± 1) = X(^,^'). 

In the following numerical simulations, we assume that the random impulses are generated 
by a Poisson process, and fix c so that all impulse marks are identical. The inter-impulse 
interval is exponentially distributed, 

W{T) = —exp(-—], (12) 



Tp \ TpJ 

where the parameter Tp is the mean impulse interval. We further simplify the calculation by 
neglecting the dependence of -R(7j, r) on r in Eq. ([9]) by replacing it with -R(7i, Tp), a normal 
distribution with fixed variance e^rp, which is equal to the average variance of the diffusion 
7j in a mean inter-impulse interval Tp. Defining G'_ = G{ip' + C,'/^, c) — G{ip' — ^'/2, c) and 
G' = G(^' + f/2,c) + G{ij' - f/2,c), the function T{^,^,^',C) can then explicitly be 
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calculated as 

peven ^ ' 

where erf is the Gauss error function. In numerical calculations, using the first several terms 
in the summation for p is sufficient. Since the error function approaches 1 (—1) very quickly 
for positive (negative) values of its argument, for a small enough value of D, the sum over 
g is to a good approximation a geometric series. 



III. NUMERICAL SIMULATIONS 

As an example of a limit-cycle oscillator, we employ the FitzHugh-Nagumo (FHN) neural 
oscillator 29| driven by common Poisson impulses and independent Gaussian-white noises 
described by the following set of equations: 

iii = e{vi + a - hui), 

3 N{t) 

^ - + Jo + a{vi, c) hit - t„) + y^r]^{t). (14) 



n=l 



Here, parameters e, a, h are fixed at e = 0.08, a = 0.7, b = 0.8, and we use the parameter 
Jo as a bifurcation parameter. The last two terms of the equation for v describe impulses 
and noises, where h(t) represents a unit impulse and a{v, c) describes f j-dependent effect 
of the impulse to the oscillator. In this example, both H and cr have only one non-zero 
component. For simplicity, we take the impulse strength c to be a constant value. When both 
terms are zero, a limit cycle exists for Jq G [0.331, 1.419], which is created by a subcritical 
Hopf bifurcation at either limits of Iq. For the simulations, we employ Jq = 0.34 and 
Jo = 0.875, which give oscillator periods of T ~ 46.792 and T ~ 36.418, respectively. We 
choose these values because the oscillator characteristics change in such a way as to show 
synchronized and desynchronized states for additive impulses, and stable 2-cluster states for 
linear multiplicative impulses. We set the mean interval between the impulses at rp = lOT. 
Results similar to the following have been obtained using Stuart-Landau and Moris-Lecar 
oscillators. However, we restrict our discussion to the FitzHugh-Nagumo model as it displays 
all of the salient features of interest. 
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In direct numerical simulations of Eq. (fT^ . we realized the Stratonovich interpretation 
by using a colored Gaussian noise generated by the Ornstein-Uhlenbeck process Tri{t) = 
~vit) + xit): where x(t) is a Gaussian white noise of unit intensity, and delivering the 
impulses as discontinuous junips of amplitude given by the Marcus approximation theorem 
of continuous physical jumps 7|, |2l|]. The correlation time r of T]{t) was set to 0.05, which is 
much shorter than the oscillator period T. In calculating the Frobenius- Perron equation (fTOl) . 
we numerically estimate X{C,,C,') and U{^) on discrete grids of dimensions between 128 to 
2048 for ^ and depending on how rapidly X(^,^') varies as a function of ^ and 
Generally, the larger the value of D, the lower the required resolution. 

We show examples of PRCs for different values of the impulse strength c obtained for the 
FHN oscillator through simulation in Fig. [H as well as the resultant transition probability 
^'). In all of the figures, we only show G [—0.5, 0.5] as X{C,, and U (^) are periodic. 
The Lyapunov exponent A is negative for the smooth PRCs, and positive for the rapidl 
fluctuating PRCs. The generic dynamical behavior of the oscillators are as follows 
When A < 0, the system settles down into a largely quiescent state once synchronization is 
achieved. The rare but sudden disintegration of a pair of oscillators is possible if there are 
regions of the PRC with positive local Lyapunov exponent, but the relative separation of a 
pair remains largely static. However for A > 0, disintegration of a pair happens routinely, 
followed by a gradual reunion, and this cycle continues ad infinitum. These occasional 
sudden, large excursions from the synchronized state is generally known as modulational 



dly 

ft 



or on-off intermittency 30, |3l| , and is a characteristic behavior of a random multiplicative 



process, of which our system is an example. 

Now let us examine the stationary distribution U{^) of the phase difference ^. We expect 
the distribution of C, to be qualitatively different between A of different sign. Figures [2] and 
Oshow the distribution of ^ for additive impulses {<j{v, c) = c, c = 0.5, —0.2, respectively) at 
various intensities of independent noise for PRCs with negative and positive A. In all figures, 
theoretical curves obtained using our Frobenius-Perron equation for the phase difference 
nicely fit the results of direct numerical simulations, which indicates that the approximations 
we have made so far are reasonable for the parameter values we use. It is readily apparent 
that if the synchronized state is stable, the synchronized peaks become taller and narrower 
as the diffusion is made smaller, while ^ far away from the stable peaks become increasingly 
rare. On the other hand, if the synchronized state is unstable, the distribution for rare ^ 
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reaches a limiting value, while only the tip of the synchronized peak increases in height and 
the width of the peak remains constant. The distributions exhibit a power-law dependence 



near ^ = 0, a characteristic of random multiplicative processes [30|, l31|, l32|, |33|]. As shown 
in Fig. m different power-law exponents are obtained by changing the impulse strength, c 
(= —0.2,0.05,0.1), where the Lyapunov exponent A determines whether the slope of the 



30 



31 



32 



power law is steeper or shallower than —1 

Figure [5] shows the same basic mechanism at work for the case with linear multiplicative 
impulses {<j{v, c) = cv, c = 0.5), which exhibits symmetric 2-cluster states. The distribution, 
which is nicely fitted by the theoretical curve, has three peaks in this case, corresponding to 
the three possible phase differences in the 2-cluster states = and ^ = ±0.5, where ^ = 
-fO.5 and ^ = —0.5 represent the same phase difference). Near each peak, the distribution 
exhibits power-law dependence as for the case of additive impulses. 



IV. COMPARISON WITH COUPLED OSCILLATORS 

We have shown that common random impulses applied to a pair of uncoupled limit- 
cycle oscillators generally produce phase coherence. Much existing work focuses on the self- 
organizing coherence brought about through coupled elements, so we would like to touch 
upon the similarities and differences between the coherence observable between coupled sys- 
tems and uncoupled systems receiving a common random input. For simplicity, we consider 
a pair of identical oscillators. 

Sufficiently weak common random input to uncoupled oscillators always tend to stabilize 
the synchronized state at zero phase difference regardless of the shape of the PRC. The 
probability density function U{C,) of the phase difference ^ always has a peak at ^ = 0, as we 
have seen in Figs. [21 [3] and [5l When the common input is stronger, the in-phase synchronized 
state ^ = can be unstable. We nevertheless observe that U{C,) has a local maximum at 
,^ = as shown in Fig. 3 for weakly unstable situations. For much stronger inputs, the PRC 
can take highly irregular forms that contain many discontinuities or with many rapid, large 
amplitude oscillations. It is then possible for U{^) to have a local minimum at = 0. 

In contrast, for oscillators with weak mutual coupling, the in-phase synchronized state 
may either be stable or unstable depending on the shape of the PRC and the interaction 
function between the oscillators. If the in-phase state is unstable, there would be no peak 
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appearing at ^ = 0; instead, a peak would be expected at some other ^ 7^ [28|, |34|. 

This illustrates the biggest difference between coherence in mutually coupled systems and 
uncoupled systems subject to common inputs. In coupled systems, it is possible to have a 
single stable phase-locked state with ^ 7^ 0, while in uncoupled systems, this is not possible. 
One possible point of confusion that arises here may be our use of the terms "stable" 
and "unstable". For uncoupled oscillators driven by common input, these terms represent 
statistical stability of the synchronized state. Even if the synchronized state induced by 
common input is slightly unstable, distribution of the phase differences can still have a 
shallow maximum at zero phase difference. The vicinity of ^ = is an attractive region even 
if the synchronized state is weakly unstable. In contrast, these terms represent deterministic 
stability for coupled systems. If it is unstable, we never observe such a maximum even if 
independent noises are added. 

If the natural frequencies of the oscillators are different, the difference in phase coherence 
behavior will be more subtle. In this case, a local extremum in f/(^) at ^ 7^ appears for 
two non-identical oscillators driven by common input, and may be a maximum or minimum 
depending on the degree of statistical stability or instability of the locked state (data not 
shown). In weak mutually coupled systems, the deterministic stability is once again depen- 
dent on the interaction function, and in addition, the magnitude of the difference of the 
natural frequencies. Furthermore, combined effects of coupling and common input, which 
may be important in practical situations, will lead to more intriguing behavior. 

V. SUMMARY 

We have found an approximate method to calculate the steady-state probability distri- 
bution of the pair-wise phase difference in an ensemble of uncoupled oscillators receiving 
random impulses. The system is essentially a random multiplicative process, and as such 
shows modulational intermittent behavior and power-law dependence of the distribution 
near its peak. Qualitative and quantitative features of the distributions have been found 
relating the results to the Lyapunov exponents that characterized the stability of clustered 



states in earlier works 
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Our treatment is conceptually a generalization of our previous result pJ3] on uncou- 
pled limit-cycle oscillators subject to common and independent infinitesimal Gaussian-white 
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noises. In that case, the common noise always stabihzes the synchronized state as long as an 
oscillator possesses a continuous phase sensitivity function. The oscillators form one or more 
synchronized clusters, depending on the degree of symmetry possessed by the system. By 
contrast, in the scenario studied in this paper, there is the further possibility that common 
impulses may destabilize the synchronized state, which can still quantitatively be analyzed 
within our theoretical framework based on the averaged Frobenius- Perron equation 391 ]. 

In this work, we considered a pair of identical oscillators subject to the same common 
impulses, and considered the diffusion in between received impulses as the effect of inde- 
pendent noises. Our method can also be applicable if the natural frequency or the PRC of 
the oscillators are slightly different. Furthermore, we can also interpret the diffusion as the 
result of inherently noisy response of an oscillator to pulsatile inputs. The consequences of a 
noisy PRC has been treated recently in the case of mutually coupled neural oscillators |28|. 



Mildly chaotic, non-mixing oscillators also show a similar noisiness to their responses. A 
noisy PRC also arises in the case of globally coupled oscillators exhibiting a collective coher- 
ent oscillation, where the response of the collective oscillation is inherently fluctuating due 
to finite-size effects, in particular near the critical point of synchronization transition 35|. 
The method developed within this paper may prove to be useful in analyzing the dynamics 
of such systems. Further results will be reported on in the near future. 

This work is supported by the Grant-in-Aid for Young Scientists (B), 19760253, 2008, 
from MEXT, Japan. 



[1] Z. F. Mainen and T. J. Sejnowski, Science 268, 1503 (1995). 

[2] M. D. Binder and R. K. Powers, J. Neurophysiol 86, 2266 (2001). 

[3] R. F. Galan, N. F. Trocme, G. B. Ermentrout, and N. N. Urban, J. Neurosci 26(14), 3646 
(2006). 

[4] R. Roy and K.S. Thornburg, Jr., Phys. Rev. Lett. 72, 2009 (1994); A. Uchida, R. McAllister, 

and R. Roy, Phys. Rev. Lett. 93, 244102 (2004). 
[5] K. Yoshida, K. Sato, A. Sugamaga, J. Sound and Vibration 290, 34 (2006). 
[6] H. Yip, S. Sano, A. Uchida and S. Yoshimori, "Multiple basins of consistency in a Mackey- 

Glass electronic circuit driven by colored noise" in Proceedings of NOLTA 2007, Vancouver, 



12 



Canada (2007). 
[7] K. Arai and H. Nakao, Phys. Rev. E 77, 036218 (2008). 

[8] Y. Zhai, 1. Z. Kiss, P. A. Tass, and J. L. Hudson, Phys. Rev. E 71, 065202(R) (2005). 

[9] H. Ukai, T. J. Kobayashi, M. Nagano, K. Masumoto, M. Sujino, T. Kondo, K. Yagita, Y. 
Shigeyoshi and H. R. Ueda, Nature Cell Biology 9, 1327 (2007). 
[10] R. Toral, C. R. Mirasso, E. Hernandez-Garci'a, and O. Piro, Chaos 11, 665 (2001). 
[11] C. Zhou and J. Kurths, Phys. Rev. Lett. 88, 230602 (2002). 
[12] K. Pakdaman, Neural Comput. 14, 781 (2002). 

[13] J. Teramac and D. Tanaka, Phys. Rev. Lett. 93, 204103 (2004); Prog. Theoret. Phys. Suppl. 
161, 360 (2006). 

[14] D. S. Goldobin and A. Pikovsky, Phys. Rev. E 71, 045201(R) (2005); Physica A 351, 126 

(2005); Phys. Rev. E 73, 061906 (2006). 
[15] K. Nagai, H. Nakao, and Y. Tsubo, Phys. Rev. E 71, 036217 (2005); H. Nakao, K. Nagai, and 

K. Arai, Prog. Theoret. Phys. Suppl. 161, 294 (2006). 
[16] H. Nakao, K. Arai, K. Nagai, Y. Tsubo, and Y. Kuramoto, Phys. Rev. E 72, 026220 (2005). 
[17] H. Nakao, K. Arai and Y. Kawamura, Phys. Rev. Lett. 98, 184101 (2007). 
[18] R. F. Galan, G. B. Ermentrout, and N. N. Urban, Phys. Rev. E 76, 056110 (2007). 
[19] A. T. Winfrcc, The Geometry of Biological Time ( Springer- Vcrlag, New York, 2001). 
[20] Y. Kuramoto, Chemical Oscillation, Waves, and Turbulence (Springer- Verlag, Tokyo, 1984) 

(republished by Dover, New York, 2003). 
[21] S. I. Marcus, IEEE Transactions on Information Theory IT-24, 164 (1978). 
[22] R. A. Gray and N. Chattipakorn, Proc. Natl. Acad. Sci 102, 4672 (2005). 
[23] R. F. Galan, G. B. Ermentrout and N. N. Urban, Phys. Rev. Lett. 94, 158101 (2005). 
[24] T. Tateno and H. P. C. Robinson, Biophysical Journal 92, 683 (2007). 

[25] F. B. Hanson, Applied Stochastic Processes and Control for Jump-Diffusions (SIAM Books, 
2007). 

[26] A. Lasota and M. C. Mackey, Probabilistic properties of deterministic systems (Cambridge 

University Press, Cambridge, 1985). 
[27] E. Ott, Chaos in Dynamical Systems (Cambridge University Press, Cambridge, 2002). 
[28] G. B. Ermentrout and D. Saunders, J. Comput. Neurosci. 20, 179 (2006). 
[29] C. Koch, Biophysics of Computation (Oxford University Press, Oxford, 1999). 

13 



[30] H. Fujisaka and T. Yamada, Prog. Theor. Phys. 69, 32 (1983); H. Fujisaka, Prog. Theor. 

Phys. 70, 1264 (1983); H. Fujisaka and T. Yamada, Prog, of Theor. Phys. 74, 918 (1985). 
[31] A. S. Pikovsky, Phys. Lett. A 165, 33 (1992). 
[32] H. Nakao, Phys. Rev. E 58, 1591 (1998). 
[33] S. Kitada, Physica A 370 539, (2006). 

[34] B. Plenty, G. Mato, D. Golomb, and D. Hansel, Neural Comput. 17, 633 (2005). 
[35] Y. Kawamura, H. Nakao, K. Arai, H. Kori, and Y. Kuramoto, Phys. Rev. Lett. 101, 024101 
(2008). 

[36] For an ensemble of uncoupled oscillators, no many-body effects due to coupling arise, and 
analyzing the phase relation between 2 oscillators is sufficient to understand the situation for 
N oscillators. 

[37] For the Ito interpretation of the impulse term, the PRC is simply given by G{4>, c) = (P{Xq{(P) + 
o-(Xo(0),c))-</> U. 

[38] The Stratonovich interpretation of Eq. ^ introduces a phase-dependent drift term that dis- 
appears upon averaging over the limit-cycle |l7], so the additive diffusion term 7^"''' may be 
taken to be zero mean. 

[39] The slope of the power-law dependence of U (.^) near the peak is always —2 for the infinitesimal 
Gaussian-white drive, while it can take a range of values in the present impulsive drive. 



14 




" -0.02 H 

o 

-0.04 





= 0.1 


c 


= 0.3 


c 


= 0.5 



0.25 



0.5 



0.75 



-0.5 



0.5 



FIG. 1: (Color online) a) The PRC G{(j)) for various values of additive impulse intensity c for 
the FHN oscillator with / = 0.34, with the PRCs of smaller amplitudes shown enlarged in the 
inset, b) The averaged phase difference transition probability for additive impulses with 

c = —0.2, D = 2.5x 10~^, corresponding to the case shown in Fig. [2l c), d) The PRCs for / = 0.875 
with multiplicative impulses, and the corresponding transition probability for c = 0.5, D = 2.5 x 
10^^, corresponding to the case shown in Fig. O The PRC of FHN gains additional symmetry 
= G(0+O.5) (as does the transition probability X{^, = X(^ib0.5, ^'±0.5)) with application 
of balanced, multiplicative noise, cr{v,c) = cv. 
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FIG. 2: (Color online) Comparison of {7(^) for the case of A < calculated using the averaged 
Frobenius- Perron equation (FPE) and measured via simulation (Sim), a) shows the global distri- 
bution in semi-log scales, and b) shows the distribution near ^ = in log-log scales for ^ > 0. The 
intensity of independent, additive noise (diffusion) is varied (D = 9 x 10^*^, 1 x 10^^, 2.5 x 10~^) 
while the intensity of the common impulses (c = 0.5) is kept constant for FHN oscillators with 
Iq = 0.875. It can be seen that lowering the independent noise narrows and increases the height of 
the peaks of the distribution near ^ = 0. Because the Lyapunov exponent remains constant, the 
slope is preserved for various diffusion strengths. 
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FIG. 3: (Color online) Comparison of for the case of A > calculated using the averaged 
Frobenius- Perron equation (FPE) and measured via simulation (Sim), a) shows the global dis- 
tribution in semi- log scales (note the y-axis range in comparison with Fig. [5] and Fig. [S]), and b) 
shows the distribution near ^ = in log-log scales. The intensity of independent, additive noise is 
varied {D = 9 x IQ-^, 1 x IQ-^^ 2.5 x 10"^) while the intensity of the common impulses (c = —0.2) 
is kept constant for FHN oscillators with Iq = 0.34. Due to the inherent instability of the ^ = 
state, the distribution of ^ reaches a limiting value as the independent, additive noise is lowered. 
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FIG. 4: (Color online) Power-law distributions of phase difference C/(^) near ^ = in log-log scales 
for the FHN oscillator with / = 0.34. The intensity of independent, additive noise is kept constant 
{D = 1 X 10^^) while the intensity of the common impulses are varied (c = —0.2, 0.05, 0.1). As the 
Lyapunov exponent of the system is changed, the slope of the power-law changes correspondingly. 
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FIG. 5: (Color online) Comparison of 2-clustered ^ distribution for the case of A < calcu- 
lated using the averaged Frobenius-Perron equation (FPE) and measured via simulation (Sim) for 
impulses with c = 0.5, FHN bifurcation parameter /q = 0.875 and independent additive noise 
(L> = 9 X 10~®, 1 X lO"'^, 2.5 X 10"^). a) shows the global distribution in semi-log scales, and b) 
shows the distribution near ^ = in log-log scales. 
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